A MILTISCALK SPLINE DERIVATIVE-BASF!) TRANSFORM FOR IMAGE FUSION AND ENHANCEMENT By IZTOK KOREN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1996 Copyright 1996 by Iztok Koren Mojim starsem ACKNOWLEDGEMENTS I would like to thank my mentors Dr. Fred J. Taylor and Dr. Andrew F. Laine for their help and support. They gave me the freedom to explore research directions that interested me most and provided me with the excellent work environment. I am grateful to Dr. Jose C. Principe, Dr. John M. M. Anderson, and Dr. Kermit Sigmon for serving on my committee and for their comments on my disser- tation. I would also like to thank all the members of the High Speed Digital Architec- ture Laboratory and the Image Processing Research Group for their friendship and help. Last but not least, I want to thank my family for their support and under- standing, and for bearing with me when I was unbearable. IV TABLE OF CONTENTS ACKNOWLEDGEMENTS iv LIST OF TABLES vii LIST OF FIGURES viii ABSTRACT ix CHAPTERS 1 INTRODUCTION 1 1 . 1 Shortcomings of Traditional Methods of Wavelet Analysis .... 1 1.2 Thesis Overview 4 1.3 Notation 5 2 SIGNAL PROCESSING USING CENTRAL B-SPLINES 8 2.1 Central B-Splines: Definition and Properties 8 2.2 B-Spline Signal Interpolations 9 2.3 B-Spline Signal Approximations 14 3 TRANSFORM IN ONE DIMENSION 17 3.1 1-D Discrete Dyadic Wavelet Transform Revisited 17 3.2 Remarks 28 4 TRANSFORMS IN TWO DIMENSIONS 30 4.1 2-D Discrete Dyadic Wavelet Transform Revisited 30 4.2 Steerable Functions 36 4.3 Steerable Dyadic Wavelet Transform 38 4.4 Multiscale Spline Derivative-Based Transform 40 5 IMPLEMENTATION ISSUES 49 5.1 Finite Impulse Response Filters 49 5.2 Infinite Impulse Response Filters 55 6 APPLICATIONS 59 6.1 Image Fusion 59 6.2 Image Enhancement 61 7 CONCLUSION 66 REFERENCES 68 BIOGRAPHICAL SKETCH 72 VI LIST OF TABLES 2.1 Transfer functions of direct B-spline filters for orders from to 9. . . 12 3.1 Impulse responses h(n), g(n), l(n), and k(n) for p = and d € {1,2,3}. 25 3.2 Impulse responses h(n), g(n), /(??), and k(n) for p=l and d € {1,2.3}. 25 3.3 Impulse responses h(n), g(n), /(n), and k(n) for p — 2 and c/ € {1,2,3}. 26 4.1 Impulse responses tx(n) for p € {0,1,2} 32 VII LIST OF FIGURES 1.1 Discrete wavelet transform of two signals translated to each other. . . 2 1.2 Discrete wavelet transform and "algorithme a trous." 3 2.1 Spline functions p (x) and r) p (x) for p € {0,1,2,3,4} 10 2.2 Fourier transforms 3 p (ui) and rj p (u>) for p € {0,1,2,3,4} 13 2.3 Splines /5 p (x) and /3 p (w) for p € {0,1,2,3,4} 16 3.1 Filter bank implementation of a 1-D discrete dyadic wavelet transform. 21 3.2 Comparison of two discrete implementations using sinc(.r) as an input. 23 3.3 Wavelets r(.r) = d \j (x) for p € {0. 1,2} and d G {1.2.3} 27 4.1 Filter bank implementation of a 2-D discrete dyadic wavelet transform. 35 4.2 Filter bank implementation of a spline derivative-based transform. . . 48 6.1 Comparison of image fusion results 62 6.2 Comparison of image enhancement results 64 6.3 The magnitude of filter coefficients for G$(uJ x ,ujy) 65 vin Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A MULTISCALE SPLINE DERIVATIVE-BASED TRANSFORM FOR IMAGE FUSION AND ENHANCEMENT By Iztok Koren December 1996 Chairman: Dr. Fred J. Taylor Cochairman: Dr. Andrew F. Laine Major Department: Electrical and Computer Engineering Analyzing images along distinct scales is advantageous in many image pro- cessing and computational vision tasks. Traditional wavelet transforms have become a popular tool for multiscale image analysis and image compression suffer from arti- facts, lack shift and rotation invariance, and may introduce aliasing and anisotropies. We propose a highly redundant representation that eliminates these artifacts. A new transform is constructed as an extension of the discrete dyadic wavelet transform with a wavelet being the first derivative of a central B-spline. We begin the construction in one dimension by providing an efficient initialization procedure for the computation of the discrete dyadic wavelet transform and by extending the transform to encompass higher order derivatives. This modified one-dimensional discrete dyadic wavelet transform will serve as a basis for deriving a steerable two-dimensional transform. Such a transform is IX advantageous in that it does not introduce any of the artifacts described above. How- ever, exact reconstruction is problematic via implementation in the spatial domain. We solve this problem by devising a spline-based approximation. The resulting algo- rithm for computing a multiscale spline derivative-based transform is implemented efficiently as a filter bank consisting of separable filters alone. For both finite and infinite impulse response filters present in the filter bank, we specify implementa- tion details of a realization that alleviates the boundary effects caused by the use of circular convolution. Finally, we use the derived transform for multiscale image fusion and enhance- ment. We show empirically that the transform does not introduce artifacts commonly reported for traditional wavelet-based implementations and provides more flexibility for different fusion and enhancement strategies by enabling orientation analysis. CHAPTER 1 INTRODUCTION Analyzing images across multiple scales and resolution has become a power- ful tool for solving compelling problems in computational vision, image processing, and pattern recognition. Wavelet theory encompasses multiscale and multiresolution representations, such as subband filtering [40], image pyramids [4], and scale space filtering [48], into a unified mathematical framework. In the area of image processing, there remain few research areas to which wavelet analysis has not been applied. For example, problems in image compression, denoising, restoration, enhancement, reg- istration, fusion, segmentation, and analysis, have all been approached with distinct kinds of wavelet processing. Though ubiquitous, wavelet analysis is not without problems of its own. Lack of translation invariance, one of the major problems of the wavelet transform [10], is in multiple dimensions accompanied with lack of rotation invariance. 1.1 Shortcomings of Traditional Methods of Wavelet Analysis A wavelet transform in its most commonly used orthogonal or biorthogonal forms is not translation and rotation-invariant. By translation-invariant transform, we mean a transform that commutes with a translation operator. Since we will deal primarily with discrete transforms in this work, we constrain the translation parameter to integer multiples of a sampling period. Lack of translation invariance of the discrete wavelet transform is illustrated in Figure 1.1. Here, we can clearly see how a translation of the input signal by one ■ 4 ill* II , ! i 4 in* 6 9 I r (a) (b) Figure 1.1: (a) Original signal and (b) signal translated one sample to the left with its discrete wavelet transform coefficients shown across dyadic scales 2 m , m G {1,2,3,4}. sample results in a completely different set of transform coefficients (orthogonal wavelet DAUB4 1 [10] was used in this experiment). Noninvariance under translations of an orthogonal and biorthogonal wavelet transform is due to lower sampling density at coarser scales. 2 A straightforward way of dealing with this problem is to construct a redundant transform by using the same sampling frequency for the input signal and all scales of the transform. A filter bank implementation of such a transform, called "algorithme a trous" [15], is based upon the fact that downsampling followed by filtering is equivalent to filtering with the upsampled filter before the downsampling, as shown in Figure 1.2. G(z) 2t H(z) 2t G(z) 2t G(z) G(z) 2t H(z) G(r) H(z) 2t (a) G(z 4 ) H(r) H(z') (b) Figure 1.2: Filter bank implementation for (a) a discrete wavelet transform and (b) "algorithme a trous" decompositions for three levels of analysis. Lack of rotation invariance is another shortcoming of traditional (i.e.. orthog- onal and biorthogonal) wavelet techniques. In defining rotation invariance, we are a bit less strict than with translation invariance. We do not require that the transform commutes with a rotation operator here. Even in the case of a simple filtering, this would limit us to circularly symmetric filters only. Our requirement for analysis is a lr rhe number in DAUB4 refers to twice the order of the wavelet (i.e.. two in this case). In practice, since analysis is performed over a finite range of scales, a discrete wavelet transform is translation-invariant by translations determined by the coarsest scale (e.g., sixteen samples for the analysis from Figure 1.1) [10]. transform that enables rotation-invariant processing. As an example of such a trans- form, let us consider filtering with the first derivative of a two-dimensional Gaussian probability density function in two directions, specifically, along x and along y-axis. By linearly combining the results of filtering in these two directions, filtering with the first derivative of a Gaussian in any direction can be computed. This fact was used by Canny [6] for edge detection. A determined edge direction rotates as an input image is rotated. After choosing the fundamental properties of the transform, one must decide upon the basis functions to be applied. For our studies, we selected basis functions that well approximated derivatives of a Gaussian, because (1) the Gaussian proba- bility density function is optimally concentrated in both time and frequency domain, and thus suitable for time-frequency analysis, (2) higher order derivatives of a Gaus- sian can be, similar to the first derivative, used for rotation-invariant processing [12], and (3) the Gaussian function generates a causal (in a sense that a coarse scale de- pends exclusively on the previous finer scale) scale space [2]. The last property makes possible scale-space "tracking" of emergent features. 1.2 Thesis Overview We wish to construct an invertible, translation and rotation-invariant two- dimensional transform that decomposes an image using a set of basis functions which are approximations to derivatives of the Gaussian probability density function. To obtain a translation-invariant representation, we will use "algorithme a trous," or more precisely, we will extend the discrete dyadic wavelet transform [28], which uses "algorithme a trous." Derivatives of a Gaussian will be approximated by wavelets which are derivatives of central B-spline functions. In Chapter 2, we review the basics of B-spline processing necessary to derive transforms in later chapters. Next, Chapter 3 examines a dyadic wavelet transform in one dimension. We present extensions of the discrete dyadic wavelet transform to higher order derivatives and to even order spline functions. To compute a discrete transform from a contin- uous one, the discrete computation must be properly initialized. We devise a new initialization procedure that is shown to be more accurate than the one suggested by Mallat and Zhong [28]. After the derivation of the transform, we point out relevant connections to scale-space filtering and reconstruction from edges. In case of the first derivative of a Gaussian, a rot at ion- in variant transform can be obtained by a tensor product extension of the one-dimensional transform to two dimensions. For second order derivatives, a tensor product approach can approximate Laplacian of a Gaussian, but cannot perform orientation analysis (due to its isotropic nature). In Chapter 4, we first examine the above two cases, and then develop a new multiscale spline derivative-based transform, which in addition enables rotation-invariant processing for derivatives of orders higher than two. Given that the developed transform is highly redundant, an efficient imple- mentation is extremely advantageous. This issue is addressed in Chapter 5. When filtering of a signal is performed by circular convolution, boundary effects may result. A standard way of dealing with this problem is by mirror-extending the input signal to the filter. We combine such an extension with symmetry/antisymmetry of the filters, to achieve savings in both computation time and memory. Finally, in Chapter 6, we present two sample applications: image fusion and enhancement. We demonstrate that the new transform does not cause artifacts com- monly associated with traditional wavelet techniques and is an attractive analytic tool for designing novel fusion and enhancement strategies. 1.3 Notation We use symbols TV, Z, and R for the sets of naturals, integers, and reals, respectively. L 2 {R) and L 2 (R ) denote the Hilbert spaces of measurable, square-integrable functions f(x) and f(x,y), respectively. The inner product of two functions f(x) £ L 2 (R) and g(x) £ L 2 (R) is given by (f(x),g(x)) = f°° f(x)g(x)dx. J — oo The norm of a function /(.r) £ L 2 (R) is defined as 11/11 = JJ~Jf(x)\*dx. The convolution of functions /(x) £ L 2 {R) and </(x) £ L 2 (R) is computed as /*</(-*•)= /" f(t)g(x-t)dt J — 'OO and the convolution of two functions f(x,y) £ L 2 {R 2 ) and g{x,y) £ L 2 (R 2 ) equals f*g(x,y)= / f(tx,ty)g(x-t x ,y-t y )dt x dt y . J—oo J—oo The Fourier transform of a function /(.r) £ L 2 {R) is defined as /(w)= / f(x)e-'»*dx, J — oo and the Fourier transform of a function f(x,y) £ L 2 {R 2 ) is equal to /(u/„w y ) = r r f(x,y)e-^ + ^dxdy. J —'OO J — oo l 2 (Z) and 1 2 {Z 2 ) stand for the spaces of square-summable discrete signals f(n) and f(n r ,riy). respectively. The --transform of a discrete signal f(n) £ 1 2 {Z) is defined as F(z)= £ f(n)z-\ n=—'Oo The convolution of discrete signals f(n) £ P(Z) and g(n) £ / 2 (Z) is equal to oo f*g(n)= ^2 f( m )9(n-m), m=— oc and the convolution of discrete signals f(n x ,n y ) £ l 2 (Z 2 ) and g{n x ,n y ) £ 1 2 {Z 2 ) is given by oo oo / * g(n x ,n y )= Yl Yl f( m *> m y)g{ n r- m x , n y - m y ) . m x =—oo m-u=—oo The Fourier transform of a discrete signal f(n) £ l 2 (Z) is equal to the transform evaluated on the unit circle CO F{u>)= £ /(»>-^ n , n= — oo and the Fourier transform of a discrete signal f(n x ,n y ) £ l 2 (Z ) is defined as FKw s )= £ £ /(n„n lr ) C -^- B - +< -» B »>. For later use, we define the following functions: 1. the unit impulse function r . . 1 for x = I I) otherwise, 2. the unit step function ( 1 for x > U{X) := for x < 0. 3. the rectangular function 4. the sine function rect(x) := < n , ? for x > i • / v sin(7r.r) smc(.r) := , and 1TX 5. the unit impulse sequence £( \ J 1 ^° r ?? = I otherwise, where x £ jR and n £ Z . CHAPTER 2 SIGNAL PROCESSING USING CENTRAL B-SPLINES In this chapter we briefly review fundamentals of spline processing needed for derivations in Chapters 3 and 4. 2.1 Central B-Splines: Definition and Properties Given real numbers — oo < x < x\ < ,r 2 < . . . < x m < x m+ i < oo, a function on the interval [.ro,.r m +i] is called a spline function of order p with the knot (i.e., grid point) sequence Xi, x-2. ... .r m . if it is (1) a polynomial of degree p or less in each interval [a?,-, x i+ i]. i = 0, 1, . . . m, and (2) continuous in its derivatives up to the order p-1 on the interval [.r ,.r m +i] (i.e., C p ~ 1 [x ,Xm+i})- Here, we will concentrate primarily on basis splines (B-splines), or more pre- cisely, central B-splines having knots at i £ Z for p odd and at i + | for p even [34]. Central B-splines of order p (with p+1 knots) are defined as Figure 2.1 shows 3 p (x) and their Fourier transforms (5 p {u)) for p £ {0, 1, 2, 3, 4}. A family of functions {/3 p (x — ???)} % forms a basis of S p , a space of order p spline functions with knots at i for p odd and at / + | for p even (?' G Z) [34. 35]. Except for p = 0, the basis functions {/? p (x — m)} are not orthogonal. Let us list some properties of functions j3 p [x) [3, 34]: 1. fip{x) are nonnegative functions with a support of length p+1, p-\-l times / * V 2. /# p (.r) = fio * /?o * • • • * A)(#) 5 where "*" denotes the convolution operator, or, equivalently, in the Fourier domain: 3. P P {x) = I ((*? + x) & P - X {x + \) + (*!-*) ^(x - §)), *-¥ = M j+ !)-M*-!)- Another interesting property of B-splines is the fact that they converge to a Gaussian as their order tends to infinity. Unser et al. [44] derived the Gaussian approximation i -4 where a v — J 2 ^-- Ratio n7 ,n , for example, is already well below 0.2%. v i^ lift* 1 1 Denoting by S\ a spline function space spanned by {fi p (x — m)} for p odd and by {/3 p (x — \ — rn)} for p even (subscript in Sj refers to the fact that spline functions have knots at integers), spline spaces form a nested sequence . . . C St, C . . . C S j C Sj C S^-i C . . . C S^-i C .... By orthogonalizing this basis functions a multiresolution analysis of L 2 (R), from which the Battle- Lemarie wavelet bases stem, can be built [10]. Here, we will not pursue constructions of orthogonal, semi- orthogonal, or biorthogonal spline-wavelets any further — reader looking for a detailed treatment of this subject may find a good starting point in books [8, 9]. 2.2 B-Spline Signal Interpolations Unser et al. developed a fast digital filtering scheme for B-spline signal pro- cessing [43]. They defined a discrete B-spline of order p and expansion factor (spacing between knots) m as b p ,m(n) := ,8 P l—j , n,m € Z. 10 p-n -3 -2 -1 I a | (b) Figure 2.1: Spline functions (a) j3 p (x) and (b) t} p (x) for p <E {0, 1,2,3,4}. II Henceforth, if the distance between knots equals one, we will write b p (n) instead of Interpolation of a discrete signal s(n) G l 2 (Z) by s p (x) G S p using central B-splines s(n) = s p {x) J2 c(i)/3 p (x-i) 5 x=n (2.1) can now be written as a convolution sum s(n) = c*b p (n). (2.2) If s(n) are samples of a function s(x) bandlimited to [— k. k] (i.e., the support of its Fourier transform S(uj) is in [— 7r,7r]), it can be shown that s p {x) — ► s(x) as p — > oo [1, 36]. In [13]. they refer to a linear operator by which B-spline coefficients c(n) can be obtained from samples s(n) as a "direct B-spline transform." Equation (2.2) therefore represents "indirect B-spline transform'' of a sequence {c(n)}. After taking the --transform of (2.2), the direct B-spline filters are found to be ^p 1 ! 77 ) = Z~ 1 {[Bp( z )]~ 1 }- Since B-spline functions j3 p (x) are compactly supported, indirect B-spline filters b p (n) are finite impulse response (FIR) filters, while direct B-spline filters ^(n) are infinite impulse response (IIR) filters. Aldroubi et al. [1] showed that IIR filters b~ l {n) are stable (i.e., the region of convergence of B~ l {z) includes the unit circle [30] ) for any order p. Note that both indirect and direct B- spline filters are symmetric, which follows from the fact that central B-splines fi p {x) are symmetric. Table 2.1 shows the --transforms of direct B-spline filters for the first ten orders. We postpone the discussion on implementation details of B-spline filters until Chapter 5. L2 Table 2.1: Transfer functions of direct B-spline filters for orders from to 9. vo 8 6 z+4+z -1 384 ~? +76J+230+76J- 1 +z~ 2 120 j 2 +26-+66+26;- 1 +.— 2 46080 '4-722; 2 + 10543j4-23548+10.543j- 1 +T22j- 2 +j- 3 5040 j 3 + 120; 2 + 1191 i+2416+l 191 z~* +120z~ 2 +z- 3 10321920 : 4 +6552j 3 +331612- 2 +24852883+4675014+2485288;- 1 +331612z- 2 +6552^- 3 + ;- 4 3(12880 : 4 +502- 3 + 14608s 2 +88234z+156190+88234z- 1 +14608;- 2 +502:- 3 + -- 4 Instead of using B-spline interpolation as given by (2.1) it is sometimes con- venient to express the interpolating function s p (x) in terms of discrete samples s(n) s p( x ) = Y, s ( l )Vpi x -i), (2.3) where r/ p (.r ) is the cardinal spline of order p. In the frequency domain, cardinal splines converge to an ideal lowpass filter with cutoff frequency tv (i.e., r] p (x) — ► sinc(.r)) as p tends to infinity [1. 45], which establishes the asymptotic equivalence with Shannon's sampling theorem [37]. Using (2.1) and (2.2) with (2.3) cardinal splines can be related to B-splines: V P (x)= E b;\i)P P (x-i). (2.4) i=—oo Cardinal splines T) p (x) and rj p {w) for p G {0, 1,2,3,4} are shown in Figure 2.2. 13 -2s 2s -2» 2s -8s -6s ^ -2s 2s 4* 6a 8n -6s -4s -2s 2s 4™ 6* 8s (a) -6s -4* -21 2s 4: 6x 8* -8s -6s -4s -2 7 2* 4s 6s 6* -8s -6s -4* -2* 2s 4s 5j -6s -6s -4s -2s 2s 4s 6s 8s -8s -6s -4s -2s 2. 4s 6s 8s (b) Figure 2.2: Fourier transforms of spline functions (a) 3 p {lj) and (b) rjJu) for p £ {0,1,2,3,4}. II 2.3 B-Spline Signal Approximations Central B-splines are also simple to use when the goal is function approxima- tion. Least-squares B-spline approximation of s(x) £ L 2 {R) is achieved by computing the orthogonal projection of this function onto S p . We have ~s(x) = s p (x) = f] d(i)/3 p (x-i) (2.5) i=—oo with d(i) = (s(x)J p {x-i)), (2.6) where °° &(*) = £ a&iCO&C* - != — OO O is the dual spline of order p [45]. Spline functions /3 p (x) and /? p (:r) satisfy the biorthog- onality condition o ((3 p (x - m), i3 p (x - n)) = 6(m - n), m, n e Z. (Note that, since both /3 p (x) and ft p (x) form a basis of 5 P , they can be interchanged in (2.5) and (2.6).) O o Figure 2.3 shows functions f3 p (x) and their Fourier transforms 3 p (uj) for p 6 {0,1,2,3,4}. An interesting alternative to the minimum Z 2 -norm (i.e., least-squares) ap- proximation of a signal is obtained by computing oblique instead of orthogonal pro- jection of the signal onto the spline function space. Unser and Aldroubi proposed an independent specification of the sampling and approximation spaces [42]: a linear op- erator maps coefficients of the input signal expansion over sampling space basis into the coefficients of the approximation space basis expansion from which the signal's projection onto the approximation space is recovered. Constraining the entire system to be linear, shift-invariant for integer translations, and consistent (i.e., the system acts as an identity operator for functions that belong to the approximation space), l.-J the obtained solution for the signal approximation is the projection of a signal onto the approximation space perpendicular to the sampling space [42]. Analogously to (2.5) and (2.6) this projection can be expressed as sr(x)= f] a{i)0 r (x-i) (2.7) i=— oo with a(i) = q£ *(s(x),p 3 (x-i)), where q^} is the convolution inverse of the cross-correlation sequence ?«(») = {P$(x - i),Pr(x)) = b s+r+1 {i). When the sampling space S s and the approximation space S r are identical (i.e., 5 = r), an orthogonal projection given by (2.5) and (2.6) results. (Note that the described oblique projection is not restricted to spline function spaces — the only requirement is that both sampling space basis and approximation space basis are Riesz bases of the corresponding function spaces [42].) Signal approximation (2.7) is particularly attractive in situations where the sampling space is given a priori (e.g., by the impulse response of the acquisition device [42]) or when such an approximation is close to the optimal least-squares solution but simpler to implement than the orthogonal projection (e.g., [47]). I(i -15 -to -5 -15 -tO -5 S 10 IS 5 10 5 tO 15 10 15 0-1 s ■ -8* -6« -4* -2i ir P-2 . _J , -At 2* 4» &> Si -<*I -2* 2i 4i 6i b 5 10 15 (a) o Figure 2.3: Spline functions (a) /3 p (x) and (b) their respective Fourier transforms 4(«)forp€ {0,1,2,3,4}. CHAPTER 3 TRANSFORM IN ONE DIMENSION In this chapter, we will augment the one-dimensional discrete dyadic wavelet transform [28] to obtain a strong foundation for derivations of two-dimensional trans- forms in Chapter 4. 3.1 1-D Discrete Dyadic Wavelet Transform Revisited A discrete wavelet transform is obtained from a continuous representation by discretizing dilation and translation parameters such that the resulting set of wavelets constitutes a frame. The dilation parameter is typically discretized by an exponential sampling with a fixed dilation step and the translation parameter by integer multiples of a fixed step [10]. Unfortunately, the resulting transform is variant under translations, a property which makes it less attractive for image analysis (cf. Section 1.1). As we have already mentioned in Chapter 1, sampling the translation param- eter with the same sampling period as the input function to the transform results in a translation-invariant, but redundant representation. The dyadic wavelet transform proposed by Mallat and Zhong [28] is one such representation. Let us begin with a brief review of properties of the dyadic wavelet transform as described in [28], but included here for completeness. The dyadic wavelet transform of a function s(x) € L 2 (R) is defined as a sequence of functions nM^)} rae z- (3-1) IT IN where W m s(x) = s * i^ m (x), and m (.T) = 2 -m ^(2 -m x) is a wavelet tp(x) expanded by a dilation parameter (or scale) 2 m . Note the use of convolution instead of an inner product. To ensure coverage of the frequency axis the requirement on the Fourier trans- form of V'm(-c) is the existence of A\ > and B\ < oo such that oo A, < £ |0(2 m cu)| 2 < B x m=— oo is satisfied almost everywhere. The constraint on the Fourier transform of the (nonunique) reconstructing function \(x) is m=— oo A function s(x) can then be completely reconstructed from its dyadic wavelet trans- form using the identity oo S ( X ) = Yl W ^ S * Xm(x), m=— oo where \ m {x) = 2- m \(2- m .r). In numerical applications, processing is performed on discrete rather than continuous functions. When the function to be transformed is in the discrete form, the scale 2 m can no longer vary over all m £ Z. Finite sampling rate prohibits the scale from being arbitrarily small, while computational resources restrict the use of an arbitrarily large scale. Let the finest scale be normalized to 1 and the coarsest scale set to 2 . The smoothing of a function s(x) G L 2 {R) is defined as S m s(x) = s*<f> m (x), where <f> m (x) = 2~ m (p(2~ m x) with m G Z, and <j>(x) is a smoothing function (i.e., its integral is equal to 1 and <f>(x) -+ as Ixl — > oo). 1!) In [28], a real smoothing function (f>(x) was selected, whose Fourier transform satisfied |^)| 2 = f)^(2 m w)x(2 m o;). (3.2) m=l In addition, it was shown that any discrete function of finite energy s(n) £ 1 2 {Z) can be written as the uniform sampling of some function smoothed at scale 1, i.e., s(n) = Sof(n), where f(x) G L 2 (R) is not unique. Thus, the discrete dyadic wavelet transform of s(n) for any coarse scale 2 M was defined as a sequence of discrete func- tions {S M f(n + a), {W m f(n + s)} m£[hM] } neZ , where s is a ir{x) dependent sampling shift. The above initialization s(n) - S f(n) is rather standard in the discrete wavelet transform computation [10]. although it yields correct results (i.e., the dis- crete wavelet transform is equal to the samples of its continuous counterpart) only when s(n) = Sqs(ti). Here, we will concentrate on wavelets which are derivatives of spline functions and this will lead us to a simple initialization procedure [46] that alleviates the above problem. For a certain choice of wavelets the discrete dyadic wavelet transform can be implemented within a fast hierarchical digital filtering scheme. Next, we shall summarize the relations between filters, wavelets, and smoothing functions. First, let us introduce a real smoothing function <p(x) such that (3.2) can be rewritten as 1 oo #(w)0M= £^(2"V) x (2 m a;), (3.3) m=0 and let us set 6{x) = /3 p (x) (i.e., we restrict ourselves to wavelets which are spline functions). ^ote that the sum index determines the range of scales of the discrete transform: using (3.2) we have ip(2u>) and \(2u) at the finest scale of the transform, while for (3.3) we get ip(u) and \(u). 20 Computing (3.3) for the finest two scales shows that 0(w) x(u) = A(«)^(«) - &(2w) 0(2w). (3.4) /5 p (2w) can be related to /i p (u;) by expressing /i p (2u>) as (cf. Proposition 1 of [46]) 1 / sin(^) \ / sin UJ p (2w) 2P+1 \™(l). and using ELo^'"^ = !i ^lf 1 e4 ^ ): 4,(2w)=(oob(|)) P+1 A(«). (3.5) (Note that a relation similar to (3.5) can be derived for integer scales provided that the dilation parameter and the order p are not both even [46].) Let F(m) be a digital filter frequency response and let F s (uj) = e JUJS F(io). If we choose 0(w) = G_ # (w)4p(«), (3.6) tp(2u)=L,(u)<p(u) t (3.7) x(u) = K.(u)<p(w), (3.8) and H(u:) = e^(cos(f)Y'\ (3 " 9) where s G {0. |} is a filter dependent sampling shift needed for g(n). /(/?). k(n), and h(n) to be FIR filters, and insert Equations (3.5)-(3.9) into (3.4), we observe the relation between frequency responses of the filters G(u)K(u) + H(u)L(u) = 1. (3.10) Similar to orthogonal and biorthogonal discrete wavelet transforms, the dis- crete dyadic wavelet transform can be implemented within a hierarchical filtering 21 scheme. To derive such a digital filtering scheme, let us assume that s(iv) from (3.1 ) is bandlimited to [— 7r,7r]. Using Shannon's sampling theorem [37] and (3.6) in the definition of the dyadic wavelet transform (3.1) with m = shows II u /oo °° °° J2 s ( i )sinc( t - i ) ]T g_ s ( m ) 3 P ( x - t - m ) dt . -oo • _, _,_ By making use of the fact that the cardinal spline functions tend to the sine function as their order r approaches infinity, and employing (2.4) we can write Wo»H * S{u>)B?{u)P r (u>)$ p {a,)G-.{u,), or, by using (3.5) and (3.9), m — 1 ^{VU m .(.r)| x= J~.SV)^- 1 ^)i? P+r+ i^)G- 5 (2 m a;)n^- s (2 ! -'). (3.11) ;=o Equation (3.11) entirely specifies the discrete dyadic wavelet transform de- composition, while the reconstruction follows from (3.4)-(3.9). Three levels of a filter bank implementation are shown in Figure 3.1. (Note that the initialization is the same as the one proposed in [46] and that except for the prefiltering and postfilter- ing, this scheme is implementing "algorithme a trous.") Noninteger shifts at scale 1 for filters with s = | are rounded to the nearest integer. — ^o» BpWco) G>) - H«o) K,(a,) Q(2<0) H(2<o) G.(4<o) HJ4<:,) ^(4(0) M4co) 1^(20.) L s (2<o) M») ^,«o); 8,(0.) Figure 3.1: Filter bank implementation of a one-dimensional discrete dyadic wavelet transform decomposition (left) and reconstruction (right) for three levels of analysis. We will now perform a simple experiment which will illustrate the difference between the implementation of the discrete dyadic wavelet transform as originally 22 proposed in [28] (i.e., without prefiltering and postfiltering) and the one from Figure 3.1. Let s(x) = sinc(.r), p — 2, and g(n) = 28 (n + 1) — 28(n) (this particular choice for p and g(n) results in the same wavelet as was used by Mallat and collaborators [27, 28]). The dyadic wavelet transform of s(x) at a scale 2 m (3.1) in the frequency domain is then C*(w) = G. s {2 m u) fa(2 m u) rect (^) . (3.12) The Fourier transform of the discrete dyadic wavelet transform of s(n) = 8(n) at a scale 2 m using spline based initialization follows from (3.11) m-l F{W m s(n)} = B;\«:)B r+3 ^)G- s (2 m u,<) Y[H_ S (2'^). (3.13) i=0 while the one using the algorithm from [28] equals m — 1 F{W m s{n)} = G-,(2 m u>) J] H..(Tu). (3.14) i=0 In Figure 3.2 a comparison of magnitudes of (3.13) and (3.14) versus (3.12) is shown: in Figure 3.2(a) magnitudes of (3.12) (solid) and (3.14) (dashed) are plotted for m € {0. 1,2.3}, while the dashed curves in 3.2(b) represent magnitudes of (3.13) with r = 5. By choosing the appropriate order r, (3.13) can approximate (3.12) in the interval [ — 7r, tt] arbitrarily good, while originally proposed (3.14) has troubles at finer scales. Mallat and Zhong [28] noticed that there was a problem with their discrete transform computation, and introduced a set of constants associated with the discrete transform coefficients at dyadic scales. They chose the values of constants such that the transform coefficient modulus maxima remained constant over all dyadic scales for a step edge input signal. In relation to Figure 3.2(a) this is equivalent to multiplying T{W m s(n)} by a distinct constant for each m. Clearly, this can improve over the situation depicted by Figure 3.2(a), but can still not compete with the described spline based initialization. 23 (a) (b) Figure 3.2: (a) Fourier transform magnitudes of the dyadic wavelet transform of s(x) = sinc(x) (solid) and the originally proposed discrete dyadic wavelet transform [28] of s(n) = S(n) (dashed), (b) Fourier transform magnitudes of the dyadic wavelet transform of s(x) (solid) and the discrete dyadic wavelet transform using quintic splines for interpolation of s(n) (dashed). 24 Next, we will choose filters in the filter bank implementation of the discrete dyadic wavelet transform. As already mentioned, we are interested in wavelets which are derivatives of spline functions. G(u>) in (3.6) is therefore the Fourier transform of the difference operator centered around —3 (cf. Property 4 in Section 2): G(w) = e*" 2; sin -J . (3.15) dmod2 where d is the order of the derivative and the sampling shift for this filter is „ Since H(uj) was already given by (3.9), the remaining two filters to be deter- mined are L(u>) and K(ui). Both of them are (as is true for <p(x) and \{x)) nonunique. If we choose L(u>) such that we can express K[u) in terms of a finite geometric series having the smallest number of elements for an arbitrary p, we get L( U ) = e-^ y. (-i)- +l [ l ^ j \ ^ cos y j (3 16) and i d+1 < / //.A \ dmod'2 / P / /, ,\ \ 2ra\ L 2 ' ,-Jujs ■ I w * where [x\ denotes the largest integer smaller than x, the sampling shift for L{u)) is the same as the one for H{u) (i.e., a= ^)^ od2 ) 1 and the sampling shift for K"(w) is the same as the one for G(u)). Note that Equations (3.9) and (3.15)-(3.17) work fine from the mathematical point of view, but in practice the reconstruction may become cumbersome when both p and d are large (the lengths of impulse responses ft(n), g(n), l(n), and k(n) are p+2, d+1, (p+l)(d-(d+l)mod2)+l, and pd+(p+l)(dmod2)+l, respectively, while for the frequency responses of the decomposition filters we observe that linip^oo \H- S {uj)\ = 8 u (u + 2nir) and Yim d ->^(2j)- d |GL,(w)| = S u (u + {2n+l)ir) with n € Z). It is also worth noting that if(w) is a lowpass filter when p is even (i.e., the reconstruction function \{x) is a wavelet only for p odd). 25 Tables 3.1, 3.2, and 3.3 list impulse responses of the four filters for p 6 {0, 1, 2} and d 6 {1, 2, 3}, while Figure 3.3 shows wavelets 0(x) = ' j*d f° r the same values of p and d. Wavelets from this family have a support of length d+p+1, regularity order p, and are either symmetric or antisymmetric. Table 3.1: Impulse responses h(n), g(n), l(n), and k(n n h(n) d= 1 d = 2 d = 3 g(n) l(n) k(n) 9(n) l(n) k(n) g(») l(n) k(n) _2 i -1 0.5 1 1 -3 -0.125 0.5 -1 0.5 -0.25 -2 0.5 -0.25 3 0.625 0.0625 1 0.5 0.25 1 0.5 -1 0.625 -0.0625 2 -0.125 for p = and d € {1,2, 3}. Table 3.2: Impulse responses h(n), g{n), /(??), and k(n) for p=l and d £ {1,2,3}. n h(n) l(n) d = 1 d = 2 i/(») k(n) 9(n) k(n) -1 0.25 1 -0.0625 1 -0.0625 0.5 -1 -0.3125 -2 -0.375 1 0.25 0.3125 1 -0.0625 2 0.0625 n J = 3 fe(n) g(n) l(n) fe(n) -3 -0.015625 -2 1 -0.09375 0.00390625 -1 0.25 -3 0.265625 0.04296875 0.5 3 0.6875 0.1015625 1 0.25 -1 0.265625 -0.1015625 2 -0.09375 -0.04296875 3 -0.015625 -0.00390625 As already discussed, wavelets with p = 2 and d=\ from a family of wavelets with p even and d=\ were used in [27, 28], whereas filters with p= 1 and d = 2 from a family of filters with p odd and d = 2 were employed by Laine and collaborators [1 1 , 22, 23]. Here described transform puts no restrictions on the choice of p or d whatsoever, and uses a better initialization procedure than the one originally proposed in [28]. 2(i Table 3.3: Impulse responses h(n), g(n), /(»), and k(n) for p — 2 and d (E {1,2,3} 7? fc(n) d = 1 d = 2 9(n) l(n) k(n) 9(n) '(») fc(n) -2 0.125 -0.015625 -0.015625 -1 0.375 1 0.125 -0.109375 1 0.125 -0.125 0.375 -1 0.375 -0.34375 -2 0.375 -0.46875 1 0.125 0.375 0.34375 1 0.375 -0.125 2 0.125 0.109375 0.125 -0.015625 3 0.015625 n d = 3 &(n) g{n) l(n) fc(n) -4 -0.001953125 0.000244140625 -3 -0.017578125 0.003662109375 -2 0.125 1 -0.0703125 0.0263671875 -1 0.375 -3 0.085937 0.0908203125 0.375 3 0.50390625 0.13037109375 1 0.125 -1 0.50390625 -0.13037109375 2 0.0859375 -0.0908203125 3 -0.0703125 -0.0263671875 4 -0.017578125 -0.003662109375 5 -0.001953125 -0.000244140625 27 -0.2 -0 6 ---71- — , - ■ ]/ ~~\ \\ | ' - / 1 ' . - <' \\ b -3 -2 -1 2 3 (a) 1 - --jr- ,__ . - 05 A \\ \l V J \\ i V 0.5 - -1 V 7 \\ // - 1 5 - -2 - Figure 3.3: (a) Wavelets (a) 0(x) = ^ M , (b) 0(ar) = ^ff^ , and (c) 0(*) d3/? g; 3 3(j) for p = (dashdotted), p=l (solid), and p = 2 (dashed). 29 the zero-crossings of W m s(x) for d = 2, a scale-space image of the representation sim- ilar to [48] can be obtained. The only differences are that curves in the scale-space plane are computed for dyadic scales only and that i) m is a close approximation of a Gaussian instead of a true Gaussian. CHAPTER 4 TRANSFORMS IN TWO DIMENSIONS In this chapter, we will extend the one-dimensional transform from Chapter 3 to two dimensions 1 and derive several two-dimensional transforms with translation and at least approximately rotation-invariant decomposition. 4.1 2-D Discrete Dyadic Wavelet Transform Revisited The dyadic wavelet transform of a function s(x,y) G L 2 (R 2 ) is defined as a set of functions [28] {Wls(x,y),Wls(x,y)} m€Z , (4.1) where W^s(x,y) = s * ^ m (x,y) for i = {1,2} and ^ m (x,y) = 2- 2m t l {2- m x,2- m y) are wavelets ^'(x,y) expanded by a dilation parameter 2 m . To ensure coverage of the frequency space there must exist an A 2 > and B 2 < oc such that A 2 < £) j:iV'*'(2 m ^,2 m u; v )| 2 < B 2 m = — oo i=\ is satisfied almost everywhere. If (nonunique) functions \ l {x.y). \ 2 (x,y) are chosen such that their Fourier transforms satisfy £ ^^(2^,2^)^(2^,2-^) = 1. l = — oo z'=l the function s(x,y) may be reconstructed from its dyadic wavelet transform by 4^y)= f; s:wi«*xi.(*,y). 7n=— oo i=l where \j n (.r,y) = 2~ 2m x i (2~ m x,2- m y). 1 For extensions to higher dimensions, please refer to [18, 19]. 30 31 However, when processing discrete functions the scale 2 m may no longer vary over all m € Z . Let the finest scale be normalized to 1 and the coarsest scale set to be 2 M . Let us, similar to [28], introduce a real smoothing function (f>(x,y) such that its Fourier transform satisfies |#c*,<*)| a = £ E^(2 m ^,,2"H) \'C2 m ^2 m ^ y ). (4.2) m =0 /=1 Here, as in one dimension, a finite energy discrete function (s(n x ,n y ) £ l 2 (Z )) can be written as the uniform sampling of some function smoothed at scale 1: s(n x ,n y ) = S f(n x ,n y ), where f(x,y) G L 2 {R 2 ) is not unique, and S m f(x,y) = /* <f> m (x,y). This led Mallat and Zhong [28] to a two-dimensional analog of the one-dimensional definition of the discrete dyadic wavelet transform: 2 { S M -if{n x +s, n y +s),{W^f(n x +s,n y +s),W^f(n x +s,ny+s)} me [o t M-i]} {nxny)&Z 2- We will use, as in Section 3.1, a spline-based initialization procedure. To implement a multidimensional discrete dyadic wavelet transform within a fast hierarchical digital filtering scheme, the wavelets were chosen to be separable products of one-dimensional functions: 1 (*,y) = ^(ar)^(y), (4.3) 4>\x,y) = V(y) </>(x), (4.4) where d(.r) and 0(.r) were chosen as described in Section 3.1 (i.e.. d(x) = /3 p (x) and ^>(u>) specified by (3.6)). From (4.3), (4.4), and (3.6), we may write frfatUy) = G- s {u! r )l3 p {u T )l] p {u;y), (4.5). j> 2 (u x ,uj y ) = G- S {u y )f3 p {u x )f3 p {uj y ), (4.6) 2 As in Section 3.1, we put the finest scale of the transform at »7) = 0. 32 where G(w) is given by (3.15) for d G {1,2}. Choosing A A X (u x ,u> y ) = K a (u x )Ti(wy)P P (u) a .)P P (u) y ), X 2 (U X ,U}y) = K a (LJy)Tl(U X )P p (u) X )0 p (u}y), (4-7) (4.8) where K(u>) and Ti(uo) are digital filter frequency responses, we may compute (4.2) for the finest two scales by X> f (2-v,2^) x\2u t ,2u y ) = \Uu x ,ui v )\ 2 - |#(2w„2w,)| 3 . (4.9) , = i Inserting the terms defined by (4.5), (4.6), (4.7), (4.8), (3.5), and (3.9) with <j>{ut x ,w v ) = ^ p {u} x )^ p (ijJy) into (4.9) results in K{uj x )G{ui x )T x {u v ) + KfaWMTrfe) + \H(u x )\ 2 \H(u y )\ 2 = 1. (4.10) Equation (4.10) represents a relation between the frequency responses of the digital filters used to implement a multidimensional discrete dyadic wavelet transform and is a multidimensional analog to (3.10). Solving (4.10) for 2i(u;) by substituting K(uj)G(u)) from (3.10) yields the closed formula [28] T^)= [ -(l + \H(^\ 2 ) (4.11) In Table 4.1 we provide the filter coefficients for Ti(uj) from (4.11) for p £ {0, 1,2}. All other filters from (4.10) were already specified in Section 3.1. Table 4.1 : Im Dulse responses t i(n) for p e n p=0 p=l p=2 -3 0.0078125 -2 0.03125 0.046875 -1 0.125 0.125 0.1171875 0.75 0.6875 0.65625 1 0.125 0.125 0.1171875 2 0.03125 0.046875 •J 0.0078125 33 As in the one-dimensional case, a two-dimensional discrete dyadic wavelet transform can be implemented as a fast hierarchical filtering scheme. To derive such an implementation, we, similar to the one-dimensional case from Section 3.1, use the definition of the two-dimensional dyadic wavelet transform (4.1) and require s(uj x ,u} y ) = for \uj x \ > 7r or \u> y \ > k. Using Shannon's sampling theorem in two dimensions [16], (4.5), (4.6), and m = 0, we get oo oo /oo /-co -Q -"-' J2 Yl s(i x ,iy)smc(t x -i x )smc(t y -i y y -CO J — CO ■ l x =-rx., l y = -,yz> X • Y, 9-s{m)(3 p (x -t x - m)f3 p {y - t y ) dt x dt y n and /CO /-oo ' — ' -'-' Y. H s{i x ,i y )smc(t x -i x )smc(ty ■° oJ -° i x =~ooi v =- o oo •Pp(x - t x ) J2 9-s(m)P P (y -t y - m) dt x dt y . m — — co oo oo 111 = — CO As in one dimension, we make use of the fact that the cardinal spline functions converge to the sine function as their order r tends to infinity, and write W^s(uj x ,u y ) ~ S(u x1 u y )B~ 1 (u} x )B~ 1 (u}y)^ p+r+1 {u} x )^ p+T+1 (u; y )G- a (ut x ) an( W$s(u} x ,u y ) ~ S(u t ,u y )B r 1 (u x )B r 1 (u}y)$ p+r+ i(u x )^ p+r+1 (u} y )G- s (u; y ), or for m £ [0, M) and discrete signal processing F{W l m s(x,y) } ~ S(u x ,ujy)B;\u; x )B^(u;y)B p+r+1 (u x y x=n x ,y=n v m—1 ■B p+r+1 (u y )G^(2 m uj x ) J] H_ s (?u: x )H_ s (?uJ y ). (4.12) t'=0 and F{W 2 m s{x,y) _ _ }~S(u x ,u y )B;\u x )B; 1 (u y )B p+r+1 (u, x y 771 — 1 ■B p+r+1 (u y )G. s (2 m i, y )l[H^C2 i ^.)H_ s C2^ y ). (4.13) i=0 Equations (4.12) and (4.13) describe the decomposition part of the filter bank implementation of a two-dimensional discrete dyadic wavelet transform. The recon- struction part can be obtained from (4.7)-(4.9) with (3.5) and (3.9). The entire filter bank implementation of the transform is shown in Figure 4.1. Except for the pre- filtering and postfiltering, we readily recognize the implementation proposed in [28]. Using the fact that a wavelet V'( x ) ^ s equal to a first (d= 1) or a second {d = 2) derivative of a spline function /? p+ d(x), (4.3) and (4.4) may be rewritten as V[x,y) = — -Q^i — , i,d€{l,2], where i) i (x.y) = (3 p+d (x)My) and d 2 (x,y) = P P (x)p p+d (y). Let us denote W m s(x,y) = (W^s(x, y), W^s(x. y)). V = (£, £), A = V 2 = (|^ + Jt-), and assume that d'(x,y, ) can be approximated by O(.v.y) for both / £ {1.2}. For c/ = 1 it then follows that W m s(x,y) = 2 m V{s*# m ){x,y). (4.14) Thus for c/ = 2 we can write 2 '£Wi nS (x,y) = 2 2m A(s*# m )(x,y). (4.15) 1=1 With ${x, y) being a Gaussian, finding local extrema of (4.14) in the direction of gradient V corresponds to the filtering stage of a Canny edge detector [6], and finding zero-crossings of (4.15) corresponds to the filtering carried out with a Marr- Hildreth edge detector (Laplacian of Gaussian) [29]. (Note that both edge detectors Bftojlgffly) B p+r+1 (m x ) B^ !(cp y ) G-s(co x ) G s (co v ) H s (co x ) H^) (a) G,(2co x ) G s (2 My ) 35 H s (2co x )H s (2o) y ) K s (co x ) T^coJ 1,(0),) K s (o) ) Ks(2co x ) 1,(2(0,) T 1 (2co x )K s (2 My ) 'H ^(2^)^(2^) H -sK) H -*(co y ) (b) B p 1 +r+l (o) x )B p 1 +r fl (cQ y ) -- B r K) B r (co y ) Figure 4.1: Filter bank implementation of a two-dimensional discrete dyadic wavelet transform (a) decomposition and (b) reconstruction for two levels of analysis. H* s (u) denotes the complex conjugate of H- s (uj). 36 involve postprocessing). Edge detection based on finding local extrema of W m s(x.y) or zero-crossings of £2, =1 W^ n s(x,y) is therefore an approximation to the Canny or the Marr-Hildreth edge detector over a range of dyadic scales. The differences stem from the fact that r)(x.y) is neither a Gaussian nor is i?*(x,y) equal to i)(x,y). 4.2 Steerable Functions When extending the transform from Chapter 3 to two dimensions, one of the first questions that conic to mind is how to deal with the fact that derivatives can be defined in any direction of the plane. 3 In case of a first derivative of a Gaussian, one can simply compute first derivatives of a Gaussian in x and y directions and then determine the derivative in any direction from these two directional derivatives [6]. For higher order derivatives of a Gaussian, Freeman and Adelson [12] showed that order+1 directional operators are needed for synthesizing the operator at any orientation. In fact, functions with the property of expressing their arbitrary rotations as linear combinations of some functions are not limited to derivatives of a Gaussian. Below, we briefly restate some of the results from [12]. A two-dimensional function is called "steerable" if its rotations generate a finite dimensional space. Rotations of a steerable function f(i\0) can therefore be expressed as /M-0 o ) = £ci(0o)/,(r,0), (4.16) where 6 denotes an arbitrary angle, {cj(0 o )} stands for a set of interpolating func- tions, {/,(r, 0)} is a set of basis functions, and r = y/x 2 + y 2 and 6 = arg(.r,y) are polar radius and angle, respectively. If /(r, 6) represents a filter kernel, the result of filtering with a rotated fil- ter f{r\6 - 6 ) can be computed simply by {c;(0 o )} weighted linear combination of outputs from basis filters {f t (r,8)}. Only the outputs from basis filters need to be 3 Second derivatives of central B-splines can be used, as we saw in Section 4.1, to approximate Laplacian of a Gaussian for approximately rotation-invariant, although not directional, processing. M precomputed and then the output from a filter rotated by any angle 9q can be found by interpolating between them. When a large number of rotations of a template filter is required, the above scheme can lead to substantial savings in both computational cost (time) and memory requirements (space). The necessary condition for a function j(i\9) to be steerable is that f(r,9) is bandlimited in its polar angle: N f(r.9) = £ a n (r)e> n6 n=-N it.17) This can be verified by inserting (4.17) into (4.16) and by assuming, for con- venience, that fi{r,9) = f(r,9 — 9,): i a n (r)e- jn6 ° =Y / c i (0o)a n {r)e- jnBi , (4.18) i=i where {#,-} is a set of user defined angles and n £ [— A", 0]. 4 Since only nonzero coefficients a n (r) are of interest, both sides of (4.18) can be divided by a n (r). This yields a matrix equation from which interpolating functions c 4 (# ) can be determined: J6o .3N60 =Jd JNt cJV2 o]N6 2 ■ J" i JNBj ci(0o) c 2 (0o) ci(9 ) (4.19) For coefficients a n = the rows corresponding to each n were removed from the matrix formulation shown in (4.19). For this system to have a solution, the angles {9;} must be chosen such that the columns of the matrix are linearly independent. In [12] they proved that the minimum number of basis functions /,(;'. 9) needed to steer /(r, 9) according to (4.16) is equal to the number of nonzero coefficients a n (r) in the Fourier series expansion (4.17). To conclude this brief description of steerability, let us only remark that func- tions which are not steerable (i.e., do not have a finite number of terms in (4.17)) 4 Note that the constraints are the same for n £ [-N, -1] and n£[l, .V], so that a subset of all possible values for n € [-N, N] can be taken. 38 can be approximated with steerable functions (a singular value decomposition was employed for approximating such functions efficiently [31]), and that the technique of expressing transformed versions of a function as linear combinations of a fixed set of basis functions is not limited to rotations (extensions to translations [39], scalings [31, 39], and general transformations [14] have been reported). 4.3 Steerable Dyadic Wavelet Transform Returning to the question from the beginning of Section 4.2, the answer seems obvious: one needs to construct a steerable analog to the one-dimensional transform from Chapter 3. Steerable transforms are nothing new — quite a few [13, 17, 38, 39] have been devised, some of them [17, 38] exactly for the computation of directional derivatives. Here, we are not interested in any directional derivatives: we want to use derivatives of central B-splines which, as the order of B-splines increases, tend to derivatives of a Gaussian. We define a steerable dyadic wavelet transform of a function s(x, y) G L 2 {R 2 ) at a scale 2 m , m 6 Z, as [21] Wi a s(x,y) = 8*tj>l(x,y), (4.20) where ^(x,y) denotes ^ m (x,y) rotated by 9 U >:\„{x.y) = 2" 2m 0(2-"V,2- m y), 0(x,y) is a steerable wavelet that can be steered with I basis functions, and t = ^j-k with i e {1,2,...,/}. Analogously to the one-dimensional case (cf. Section 3.1) we require the two- dimensional Fourier plane to be covered by the dyadic dilations of ^ i {2 m u x ,T n u y ): there must exist ,4 3 > and B 3 < oo such that 4» < £ El^'(2 m ^,2 m u;,)| 2 < i?3 (4.21) m=— oo i=l is satisfied almost everywhere. 39 If (nonunique) reconstructing functions \^(:r,y) are chosen such that their Fourier transforms satisfy f) E ( / > '(2 m u; x ,2 m ^)x ! (2 m ^,2 m u; y ) = l, (4.22) m=— oo t=l the function s(x, y) may be reconstructed from its steerable dyadic wavelet transform by oo / »{*,V)= E EwL**x\n(*,v)< (4-23) m=— oo i=l where \ l m {x,y) denotes \ m {x,y) rotated by 6, and \ m (x.y) — 2~ 2m \(2 _m .r, 2~ m y). To derive an algorithm for the fast computation of the transform, we, similar to (3.3), introduce two smoothing functions such that 00 I i{L> x ,u> y )<p{u x ,u> y ) = £ ^^'(2 m ^,2 m a Jy )i'(2 m u; x ,2 m u; y ). (4.24) We choose wavelets that are steerable analogs to the one-dimensional wavelets from Section 3.1: 5 il>(u r ,ue) = (ju} r cos(u e )) d (— s^J , (4-25) where u T = «/u;J + ^ and u^ = arg(cj :r ,u; y ). These wavelets can be steered with d+\ basis functions (i.e., / in (4.16) is equal to d+1). Choosing 4>(2u r ) - H st (u r )<j)(u) r ), (4.26) ^ % {u} r ,u 6 ) = G st {uj r ,u: 9 - 9i) 6{u: r ). (4.27) ip(2u T ) = L st (u> r )ip(u r ), (4.28) and X % (w r ,u e ) = K at (L> r ,u> g - Oi)<p{u) r ), (4.29) This choice can be viewed as an extension of the transform presented in [20, 21, 24]. 10 and using (4.26) through (4.29) with (4.24) computed for the finest two scales, we obtain ^G a t(u; r ,u>4, - O^Kstiur,^ - 0i) + H 3t (ur)L s t{u T ) = 1. (4.30) Setting <£(u> r ) = /? p (uv), and employing (3.5) and (4.25) with (4.26) and (4.27), we find that H«(u r ) = H-fa) (4.31) and G at {u r ,u e ) = (co S (u e )) d G- 3 {u> r ), (4.32) where H(u>) and G(u) are specified by (3.9) and (3.15), respectively. By inserting (4.31) and (4.32) into (4.30), the missing two filters can be chosen as L«(u r ) = L a (u r ) (4.33) and K at (u r ,U0) = — (cos{u e )) d K s {u r ), (4.34) where L{^) and K(u)) are given by (3.16) and (3.17). respectively, and Cd = Ef = i(cos(u;-^)) 2d . 4.4 Multiscale Spline Derivative-Based Transform Let us pause here for a brief assessment of the two-dimensional steerable trans- form derived so far. We have chosen steerable wavelets (4.25) which are equal to d-th order derivatives of circularly symmetric spline functions in the direction of x-axis (note that knots for these splines are circles) and we have laid a foundation for fil- ter bank implementations in (4.30). An obvious shortcoming of this scheme is the fact that none of the filter kernels obtained from (4.31) through ( 1.31) is compactly supported on the rectangular grid. For implementations using digital filters, one is therefore forced to approximate these frequency responses, and by doing so. (4.30) II may not hold anymore. Filters in filter bank implementations of steerable pyramids described in [17, 38, 39], for example, were designed by using various techniques for approximating desired frequency responses. None of the reported filter banks achieved perfect reconstruction and all filter kernels were nonseparable. Here, we will take a different approach. We will approximate the wavelets in (4.25) in a way that will lead to x-y separable filters in a perfect reconstruction filter bank imple- ment at ion of the transform such that the quality of approximation will improve by increasing the order of B-splines. Let us approximate wavelets from (4.25) with , v d (3 p+ d(x) nx,y) = — -jp — i3 P +d(y)- (4.35) Based on the fact that B-splines tend to a Gaussian as their order increases, it is easy to see that both wavelets (4.25) and (4.35) converge to the same functions (i.e., d-ih order derivatives of the normalized Gaussian in the direction of .r-axis) as p — * 'CO. In order to steer wavelets tj){x,y) given by (4.35) (note that steering will be only approximate, since these wavelets are not steerable), we need to find basis functions that will approximately steer ij)(x,y). Computing rotations, as we did in (4.20), is not practical here, because arbitrary rotations of (4.35) cannot be expressed exactly in terms of central B-spline functions from Chapter 2. Instead, we take advantage of the property of circularly symmetric functions that rotations of their directional derivatives are equal to directional derivatives in rotated directions: d d Q c (x,y)\ d d g c (x,y) Tin dn d dni "0 where 1Z 6q stands for rotation by 6 , dec £' v) = nVg c (x,y), g c (x,y) is a circularly symmetric function, and ng denotes vector n = (cos 0, sin 6) rotated by 6 . Let us choose g(x,y) = 8 p+d (x)/3 p+d {y), 12 which is approximately circularly symmetric function for higher order splines. A rotation of %l>(x,y) = f^J y from (4.35) by 6q can therefore be approximated by where n = (cos $o, sin #o) = (n x ,n y ). Note that in case of Gaussian, which is both x-y separable and circularly symmetric, (4.36) becomes exact (e.g., for g(x,y) = e~^ x + y ~>, o = —0, and d = {2,4}, we obtain, up to a scaling factor, x-y separable basis set for the second and fourth derivative of Gaussian from Tables III and VII of [12]). Having found a set of basis functions (4.36) that approximately steer wavelets (4.35), we want to construct a transform such that Equations (4.20) through (4.24) will be valid (superscript i must be viewed now as an index, rather than rotation by 0i). In frequency domain, we can express basis functions from (4.36) as 4> i+1 (u x ,u; y ) = G±?(u x )GL,(Ljy)P P+i {u x )P P+i - i (u y ), £ = 0,1 d, (4.37) where G d (u) is given by (3.15) and G°(w) = 1. Choosing appropriate x % (u x ,u> y ) to obtain a relation needed for the filter bank implementation of the transform is more complicated than in one dimension. Since we could not find a general solution for arbitrary d, we solve each case separately. Below, we present solutions for the first four orders. When d < 2, we impose ^(w^u,,) = <f>(u> x ,u}y) = /3 p (u> x )fl p (u}y), a constraint analogous to the one from Section 3.1. For d = 1 , we write similar to [28] x\u>*,Uy) = K]{u x )T l {uj y )i3 p (uj x )fip-\{uy). (4.38) \ 2 (u x ,^ y ) = r 1 (w x )j£(<«;J&_ 1 (w x )&(a; y ), (4.39) where K d (u>) and Ti(u) are given by (3.17) and (4.11), respectively. 13 Computing (4.24) for the finest two scales and inserting (3.5). (4.37), (4.38), and (4.39) yields a relation between frequency responses G\u x )K\u> x )Ti{u>y) + TiK^Vv^Vv) + \H(u t )H{uy)\ 2 = 1. For d = 2, we choose \\u; r ,uj y ) = K*(u x )T 2 (uJy)Pp(u x )P P -2(u y ), (4.40) X 2 {u x ,u v ) = K](u x )K}((jj y )P P - 1 (u g )P P -.i(u y ), (4.41) X 3 (w*,w y ) = T 2 (u; ;r )A' s 2 (u;j / )^ p _2(u; ;c )^ p (LJ y ), (4.42) w here r a (w) = |#(u,)| 2 . (4.43) Using (3.5), (4.37). and (4.40) through (4.42) with (4.24) results in G 2 (u x )K 2 (u x )T 2 (u y ) + G\u; x )K\u x )G l {w y )K\u y ) + T 2 (u x )G 2 (u y )K 2 (u y )+ +\H{u x )H(u> y )\ 2 = l. For orders d > 2, we require <f>(u> x ,u y ) = J3 p (u} x )P P {u> y ) and ip(u x ,u> y ) = (f>(u} x )<p(uj y ), where fi{u>) is specified by (3.7) and (3.16). For d = 3. we choose reconstructing functions ^(u^Uy) = K*(w x )<p(u x )<p- 3 (L,j y ), (4.44) X 2 (w„w y ) = -A' s 2 (^)A' s 1 (cj y )V3(u; r )V3(^)^-i(^ x )^-2(^), (4-45) fK,w s ) = -A s 1 (^)A 2 (u; y )V3(a J;r )l / 3K)^-2(^)^-iK). (4.46) X*(ws,w y ) = K*(u y )<p- 3 (u x )<f(u) y ), ( 1.47) where 7 3 (w) = -^(1 - \H(lo)\ 2 ), (4.48) and £_,-(a>) € L l (R) denotes a function such that <p(uj) = /3,_i(u;)<£_ t -(a;), i € iV. II Employing (4.37), (3.5), (3.7), and (4.44) through (4.47) with (4.24) yields a relation -G\u, x )K\u; x )V 3 (u; x )G\u: y )K 2 (u y )V 3 (u Jy ) + G 3 MK 3 M+ +H(u x )L(u> x )H(u y )L(u y ) = 1, where L(*c) is specified by (3.16). For d = 4, our choices are X l (e>z,w y ) = Kt(oj x )T 2 {uj y )^(u) x )^- 4 M, (4.49) X 2 (wx,u y ) = K 3 (u x )Kl(u y )<p-i(u x )<p- 3 (u> y ), (4.50) X 3 (w*» w i») = - K Uvx)Kl(uy)V 4 (u x )V4(u y )<p-2(ux)<P-2(vv), (4.51) X 4 (u x ,u y ) = A' s 1 (^ x )A' s 3 (u.' y )^_ 3 (u; J .)^_i(^ !/ ), (4.52) X B (w t ,Uy) = T 2 (u> x )K 4 (u> y )<p- 4 (u x )ip(u v ), (4.53) where V 4 (u x ) = 1 - |tf(u;)| 2 . (4.54) Using the above functions (4.49) through (4.53), (4.37), (3.5), and (3.7) in (4.24) computed for the finest two scales gives G 4 (u x )K 4 (u x )T 2 (u y ) + G 3 (u; a .)A' 3 (-.v)G 1 (u; J ,)A' 1 (u; J ,)- -G 2 {u x )K 2 (u x )V 4 (u; x )G 2 {u y )K 2 {u; y )V 4 (u y ) + G\u^ +T 2 {uj x )G 4 {u; y )K 4 {uj y ) + H(u x )L{u x )H(u y )L(u> y ) = 1. Here, we have even more freedom for choosing the reconstructing functions than in one dimension. The above functions for d = {2,3,4} were found by trying to imitate the one-dimensional transform from Chapter 3 as much as possible. All decomposition filters G d (m) were first paired with corresponding reconstruction filters 45 K d (uj), and then all other potential digital filters were specified as polynomials in H- a (u). We inserted thus specified filters into a relation between their frequency responses and solved for the unknown polynomial coefficients. Since we allowed more filters with higher-degree polynomials than expected in the solution, the resulting system of linear equations was underdetermined. This allowed enough freedom for removal of undesired digital filters and for balance between degrees of polynomials. The described procedure for determination of reconstructing functions and filters involves quite a lot of heuristics to obtain the appropriate solution from the underdetermined linear system. Unfortunately, we are not aware of any systematic way (aside from numerical optimization, which may be pretty cumbersome) to obtain solutions comparable to the ones above. Next, we will derive a filter bank implementation of the transform. As in Section 4.1, we assume a bandlimited input signal: s(uj x ,u} y ) = for |o; z | > it or |a> y | > 7r. Using Shannon's sampling theorem in two dimensions [16] with (4.20) and basis functions from (4.37), we can write /OO /"CO ^ -~ J2 J2 8 (**i »v)mhc(** - i x )ainc(t y - i y )- -OO J — OO ■ l x = — OO ly = — OO OO OO Y. g d -7( m ?){3p+i(x -t x -m x ) Y. g- s {m y )f3 p+d -i{y -ty-m y )dt x dt y , m x = — oo m,u=—oo where i = 0, 1, ...,(/ as in (4.37). Again, we approximate sine functions with r-order cardinal splines, then use (2.4), and get x — n x ,y — fly m—\ ■B p+r+d . i+1 (u; y )G d - i (2 m u; x )GU2 m u y ) J] H p _ + a i (2 n u x )H^ d - i (2 n u; y ). (4.55) 71=0 Using (4.55) with an approximation 5 p+r+t+1 (w) ~ B p+r (uj)Bi(uj), we can ob- tain a filter bank implementation of the transform decomposition. The reconstruction part follows from (4.24), (4.37), and reconstructing functions for distinct values of d. 46 Figure 4.2 shows filter bank implementations of a multiscale spline derivative-based transform for d — {1,2,3,4}. For d = 1, we recognize (except for the prefiltering and postnltering) the filter bank implementation of a two-dimensional discrete dyadic wavelet transform from [28]. For d = 2, however, our transform differs from the filter bank presented in [22] (i.e., the corresponding transform described in Section 4.1): second derivative is computed only in the directions of x and y-axis in [18, 22], which is not enough for steering. Although not particularly appropriate for orientation anal- ysis, such a scheme may still, as we have seen in Section 4.1, efficiently approximate Laplacian of Gaussian across dyadic scales. A transform similar to the one described in this section, was presented in [17, 38, 39]. Their filter bank implementation is less redundant (downsampling is used on the lowpass branch, while simultaneously keeping aliasing negligible by using a filter with insignificant amount of energy for \uj r \ > |) and uses reconstruction filters with same magnitude frequency responses as the decomposition ones — a possible advantage for certain applications. They have, on the other hand, little control over the function from which derivatives are computed (to obtain a d-th derivative, they multiply a circularly symmetric filter by (-jcos6) d with all filters being obtained by recursive minimization of a weighted combination of constraints), filter bank does not have perfect reconstruction, and none of the filters is x-y separable. IT B> x ) B> y ) -JJ^hK) B p+r+1 (co y ) (a) B p+r+1 (ro x ) B p+r+1 (co y ) B r (o) x ) B r (co y ) (b) G 1 s (2 m co x ) G 1 s (2 m co v ) , sV - w y , H s (2 m (o x )H,(2V) (c) G 2 (2 m co x ) Qs(2 m a) x ) G 1 s (2 m co y ) G 2 (2 m co Y ) H,(2 m co x )H,(2V) (e) Kj2 m a> x ) T,(2 m co v ) "1 \ TT-l/^m 1,(2©,) K;(2 m co v ) L s (2 m o) x ) L s (2 m co y ) (d) K 2 (2 m (o x )T 2 (2 m co v ) •'x/ ^V* - "'y/ K;(2 m co x )K;(2 m o) y ) T 2 (2 m co x ) K 2 (2 m co v ) L s (2 m co x ) L s (2 m co y ) (f) Figure 4.2: Filter bank implementation of a multiscale spline derivative-based trans- form for m G [0, M - 1]: (a) Prefiltering, (b) postfiltering, (c) decomposition and (d) reconstruction modules for d = 1, and (e) decomposition and (f) reconstruction modules for d = 2. is G 3 s (2 m co x ) B^V) G 2 (2 m (o x )G.;(2 m (o y ) Gl(2 m co x ) G 2 (2 m co Y ) B 2 (2 m co x )Q 3 s (2 m (o y )- H s (2 m (o x )H,(2 m co y ) (g) K^coJ EfcV) -K 2 (2 m (o x )V 3 (2 m co x )K;(2 m co,)V3(2 m (o y ) K;(2 m (o x )V 3 (2 m (o x )K 2 (2 ni co y )V3(2 m co y ) ^(2 m co x )K 3 (2 m co y ) L s (2 m co x ) L S (2V) (h) G 4 (2 m (0 x ) B,(2 m co y ) G^X^coJB^V) G 2 (2 m co x ) G 2 (2 m (o v ) G^coJ B 2 (2 m a) x ) C&V) B 3 (2 m (o x ) G 4 (2 m co y ) H :S (2 m o) x )H s (23) (i) K 4 (2 m (o x )B 3 , (2 m (o y )T 2 (2V) K 3 (2X)B,\2V)K:(2 m co y ) K 2 (2 m (o x )V 4 (2 m (o x )K 2 (2 m (o y )V 4 (2 m (o y ) R 2 , (2 m a) x )K;(2 m co x )K 3 (2 m (o y ) B 3 1 (2 m (o x )T 2 (2 m (o x )K 4 (2 m (0 y ) 1 LsPX) L s (2 m co y ) 0) Figure 4.2: Continued: (g) Decomposition and (h) reconstruction modules for d = 3, and (i) decomposition and (j) reconstruction modules for d = 4. Decomposition modules are recursively applied at the locations of the filled circles. CHAPTER 5 IMPLEMENTATION ISSUES Except for the steerable dyadic wavelet transform, all transforms presented in Chapters 3 and 4 can be implemented as filter banks consisting of one-dimensional fil- ters only. In this chapter, we show how to take advantage of symmetry/antisymmetry of filters when combined with a mirror extended input signal. 5.1 Finite Impulse Response Filters Since all two-dimensional filters used in the filter bank implementations of the transforms from Chapter 4 are x-y separable, we will begin this section with a detailed description of FIR filter implementations for the one-dimensional discrete dyadic wavelet transform implementation from Figure 3.1. The extension to two dimensions will then be straightforward. Let us refer to filters applied at scale 2 m as filters at level m+1, and let filters at level 1 (Equations (3.9), (3.15) through (3.17), (4.11), (4.43), (4.48), and (4.54)) be called "original filters,'' to distinguish them from their upsampled versions. As an input to the filter bank from Figure 3.1, we consider a real signal s(n) G 1 2 {Z). n G [0, N - 1]. Depending on the length of each filter impulse response, filtering an input signal may be computed either by multiplying the discrete Fourier transforms of the two sequences or by circularly convolving s(n) with a filter's impulse response. 1 Of course, such a periodically extended signal may change abruptly at the boundaries causing artifacts. A common remedy for such a problem is realized by constructing x As is customary in image processing, we use circular, rather than linear, convolution. 1!) 50 a mirror extended signal [18] s (n) = l S{ - n - l) if « € [-iV, -1] mel j \ s(n) line [0,7V -1], ( ° j where we chose the signal s me (n) to be supported in [— TV, N — 1]. It will become evident shortly, that mirror extension is particularly elegant in conjunction with symmetric/antisymmetric filters. Let us first classify symmetric/antisymmetric real even-length signals into four types [30]: Type I f(n) = /(-»), Type II /(n) = f(-n - 1), Type III /(n) = -/(-«), Type IV /( n ) = -/(-n-l), where n € [-A r , A' - 1]. Note that for Type I signals the values at /(0) and f(—N) are unique, and that for Type III signals the values at /(0) and f( — N) are equal to zero. Using properties of the Fourier transform it is easy to show that the convolu- tion of symmetric/antisymmetric real signals results in a symmetric/antisymmetric real signal. If a symmetric/antisymmetric real signal has an even length, then there always exists an integer shift such that the shifted signal belongs to one of the above types. Now, we are ready to examine the filter bank implementation of a one-dimensional discrete dyadic wavelet transform from Figure 3.1 with filters given by (3.9) and (3.15) through (3.17) driven by a mirrored signal s me (n) at the input. Let the number of levels M be restricted by M < l + log 2 / ~* (5.2) ■'-'max -1- 51 where L max is the length of the longest original FIR filter impulse response. Each FIR filter block in the filter bank consists of a filter and a circular shift operator. Equation (5.2) guarantees that the length of the filter impulse response does not exceed the length of the signal at any block. Since our input signal s me (n) is of Type II and noninteger shifts at level 1 are rounded to the nearest integer, it follows that a processed signal at any point in the filter bank belongs to one of the types defined above. This means that filtering a signal of length 2 A* can be reduced to filtering a signal of approximately one half of its length. (For Types I and III, A + 1 samples are needed. However, for Type III one needs to store only A — 1 values because zero values are always present at the zeroth and (— JV)-th sample position). Implementation is particularly simple for FIR filters designed with d even and p odd. Filters are of Type I in this case, so the signal at any point of the filter bank will be of Type II. An FIR filter block from the filter bank shown in Figure 3.1 can therefore be implemented by F s , m u(n) = f(0)u II (n) + Y,f(i)[u II (n-2 m i) + u JI (n + 2 m t)l n € [0, A r -1], (5.3) where u//(n) = < u(-n-l) ifne[-f,-l] u(n) ifn€[0,AT-l] (5.1) u(2N-n-l) ifneJA,^-], u(n) is an input signal to a block, /(??) is an impulse response of some original filter. L is the length of the filter, and A* is the length of an input signal s(n) to the filter bank. Implementation of filters b p (n) used for prefiltering and postfiltering represents a special case of (5.3) with m = 0. A filter bank with the above implementation of blocks and signal s(n) at the input yields equivalent results as circular convolution for s me (n) as defined by (5.1). In addition to requiring one half the amount of memory, the computational savings 52 over a circular convolution implementation of blocks are, depending on the original filter length, three to four times fewer multiplications and one half as many additions. A similar approach can be used for other filters. However, things get slightly more complicated in this case, because the filters are not of the same type and the signal components within the filter bank are of distinct types. As a consequence, an implementation of blocks that use distinct original filters may not be the same, and the implementation of blocks at level 1 may differ from the implementation of blocks at other levels of analysis. The decomposition blocks at level 1 can be implemented by 2 1 G- a , u(n) = J2 g{i)[u n {n - i - 1) - u n (n + i)], n G [1, JV - 1], for d odd, (5.3) for d even, H- s ,ou{n) = Y, h{i)[un{n-i-l) + u n (n + i)], n G [O.A], t'=0 for p even, and (5.3) for p odd, where «//(/) is defined by (5.4), g(n) and h(n) are impulse responses of the filters computed from (3.15) and (3.9), respectively, and L is the length of the corresponding impulse response. The output from a block G- 3 (u>) at level 1 is of Type III for d odd and of Type II for d even, while the output from H_ a (u) at the same level is of Type I for p even and of Type II for p odd. The decomposition blocks at subsequent levels m G [l,Af — 1] can be imple- mented by G- a<m u(n) = £ g(i)[ui(n - 2 m (i + s)) - Ul (n + 2 m (i + s))}, n G [l,iV- 1], ;=0 for d odd and p even. 2 l G- a , m u(n) = Y, 9(i)[un(n - 2 m {i + s)) - u n (n + 2 m (i + s))}, n G [0, N - 1], f=0 53 L-\ 2 for d and p odd. F_ s>m «(n) = /(O)ur(n) + £ /COM" " 2™*') + "/(« + 2 "'0]- n € [0, AT], (5.5) 1=1 with f(n) = g(n) for d and p even, H.,, m u(n) = J2 h(i)[ Ul (n - 2 m (i + a)) + u/(n + 2 TO (* + a))], n € [0, AT], (5.6) i=0 for p even, and (5.3) for p odd, where »(-??) u/(n) == { u{n) ifn€[-f,-l] if n € [0, N\ (5.7 u(2N-n) \ine[N + l,m The outputs from blocks G- S (2 m u:) are of Type III for d odd and p even, of Type IV for d and p odd, and of Type I for d and p even, whereas the outputs from H- S (2 m u)) are of Type I for p even and of Type II for p odd. Next, the reconstruction blocks at level 1 can be implemented by K s .ou{n) = J2 Hi)[uin(n -?' + !)- u ni (n + i)], n € [0, N - 1], i=i for d odd, (5.3) for r/ even, L 2 L,, u(n) = £ '(»)[«/(» - »' + 1) + "/(» + *)]i n G [0, A r - 1], ! = 1 for p even, and (5.3) for p odd, where uin(n) = < -u(-n) ifne[-f. -11 if n = u(n) if n e [1,N -1] \(n = N [5.8) ( -u(2tf-n) ifne [W+l,^], «/(/?) is as defined by (5.7) and k(n) is an impulse response of the filter from (3.17). Note that both outputs from blocks K a (u) and L a (u) are of Type II. The reconstruction blocks at subsequent levels can be implemented by i-l 2 * K a<m u(n) = £ H* + !)[«///(» - 2 m (i + s)) - U//7 ( n + 2 m {i + .s))], n <= [0, AH, !=0 54 for d odd and p even, (5.5) with f(n) = k{n) for d and p even, A' s , m »(n) = £ fc(i + l)[u/v(n - 2 m (i + s))- u IV (n + 2 m (i + s))], n € [0. N - 1], i=0 for d and p odd. L sjn u(n) = #_ s , m u(n), for p even, and (5.3) for p odd, where «///(/) is given by (5.8), ' -«(-„- 1) ifne[-f.-l] u /v (n) = I n(n) if n € [O.N - 1] , -«(2JV-n-l) ifnG[A',f-], and H- S _ m u(n) is specified by (5.6). We observe that the outputs from blocks K s ('2 m uj) and L s (2 m u;), m G [1. M — 1], are of Type I for p even, and of Type II for p odd. When we compare the above implementation of blocks to circular convolution driven by a mirrored signal s me (n) at the input, we observe that approximately twofold less memory space, three to four times fewer multiplications and one half as many additions are required. (For Type I signals an additional sample has to be saved because two values are without a pair). Until now. we have talked only about the one-dimensional case, whose filter bank implementation is depicted in Figure 4.2. Two-dimensional transform filter bank implementations (Figures 4.1 and 4.2) are not only comprised of x-y separable filters solely, but also use all the filters from Section 3.1. Everything presented in this section so far is therefore directly applicable to the two-dimensional case. Filters which have nol been t reated yet (i.e.. /,(//). /_>(;/). r : ,,( n ). r,( n ). and filters b p (n ) from the decomposition modules of Figure 4.2) can all be realized by (5.3) for p odd or m = and (5.5) otherwise (/(??■) denotes an impulse response of any of the above mentioned zero-phase filters in this case). 2 In case of filters i>3(n) and 1^(71), a slightly more efficient implementation can be obtained by precomputing new filters k * i>3<n) and k * 1)4(71), and then implementing them by K s m u(n), me [0,M- 1]. I 1 u- implementation presented in this section performs all operations in the spatial domain, however, one could also implement the structures shown in Figures 3.1, 4.1, and 4.2 with an input signal s me (n) (Equation (5.1)) in the frequency domain. For short filter impulse responses, such as those given in Tables 3.1. 3.2 and 3.3, the spatial implementation described in this section is certainly more efficient. For long filter impulse responses, however, filtering is faster if implemented in the frequency domain. Additional details on alternative implementation strategies can be found in [33]. 5.2 Infinite Impulse Response Filters Implementation of IIR filters b' 1 ^) which were introduced in Section 2.2 is a bit more involved than the one encountered in the previous section. Fortunately, the number of different cases is much smaller here: possible input to b^ln) in filter banks from Figures 3.1, 4.1, and 4.2 is either of Type II or of Type I. 3 We will use ideas and a few results from [43]. Let us first take a closer look at the system function B^iz). This function can be written as a cascade of terms E ^ - j^h^ - ,i- a --.) ( i- a -, - ^ which can be expressed in a parallel form as -a / 1 1 -a 2 Vl -or 1 + Y^az E ^ = ^-i 7— nr + i-^rr-0- <wo) where a and - are poles of the causal and the anticausal filter, respectively. The impulse response of this term can be written as 1 — cr 3 Note that symmetry types in this section slightly differ from those defined in Section 5.1: here, mirror extended signals are periodically repeated, so that they stretch from — oo to oc. .-,(, We choose to implement E{z) in a cascade form and therefore extract the difference equations from (5.9): c + (n) = u(n) + ac + (n - 1) n = 1,2 V-l. (5.11) and c(n) = a(c(n + l)-c + (n)) n = N-2,N-3,. . .,0, (5.12) where u(n) denotes the input to the block, c + (n) is the output from the causal part, and c(n) stands for the output from the block. To solve (5.11) and (5.12) we need boundary conditions c + (0) and c(N — l). Let us begin with filters fcr 1 (n) in filter bank implementations from Figures 3.1, 4.1, 4.2(a). 4.2(b). and Figures 4.2(h) and 4.2(j) with m = 0. We derive ■ ' v_1 a' +l + o c+(0) = £ a-u„,(0 = «(0) + £ - T f^ :'=— oo i=0 and, using parallel form (5.10) io v27V -</(/ u(0) + 5> <+1 ti(«), (5.13) i'=0 c(N-r —a 1 -a' 2 — a 1 -a 2 ( c + (iV-i)+]T -i n -V-, + Q .v+i +i I ,2V u(iY i=0 N-l [c + (N-l)+ Yl ^-^(i)), i=.V — 1 — i'o (5.14) w here u 1Ip (n) »//(") : «//(;? mod (27V)) if n>0 »//(-(n + l) mod (2JV)) if n < 0, u(n) if n e [0, A r - 1] u(2JV-n- 1) if/? € [.V.2.V- 1], iV is the length of an input signal to the filter bank, and ?o < N — l is selected such that a la falls below a predefined precision threshold. For IIR filters from Figures 4.2(h) and 4.2(j) with m G [1. M - 1] and p odd, we get o c+ (°) = £ [a-**u IIp (i t=— oo oo N-l Y.Y. n=0 i=0 OO iV— 1 .. r -, N .„. W--V ^^ I r »+l+2Nn 1 2.Y(n+l)-i | = "(0) + EE [o^ 5 ^] + a-^r— u (i) (5.15) 57 and N-\ ^(A f -i) = T^ T (^(-v-i) + EL{h 2?n + n—0 i'=0 Q .Y+l + i + 2.Yn }«(0, (5-16) where [«1u - n x if .r € Z otherwise. If N is a power of two and 2 m_1 < JV ((5.2) guarantees that the latter condition is always true) (5.15) becomes c+(o) = «(o)+ Y. i+1 a 2 m + 2-V-i Q 2" ;=o W hile (5.16) can be written as c(X-l) = -^(c + (N-l)+Y — o ^~^ 1 — Q 2™ N-i Q 2 m a «(0, + ■ v tH' Q 2 m 1 - Q2 m -»(/)• Finally, the boundary conditions for filters ^(rc) from Figures 4.2(h) and 4.2(j) with m E [1,M — 1] and p even are c+(0)= £ [a-^u/pCi) i=— oo V^ f r 2ATn 1 ^ j [a -5™-] i/(0)+ a^ 7- + n=0 N-l + 2JVn 2A r |7i + l)- Q 2"> u(N)+ u(i)\ (5.17) and where —a c(7V) = - -(2c + (AT)-«(JV)), 1 -a uj(n) = «/p(«) = »/(H mod (2A*)). ' u(n) if n € [0. A"] u(2JV-n) if n e [N + 1,2N-1], and c~(N) = c + {N) with c~(n) denoting response of the anticausal filter from (5.10) was used. oS Again, if N is a power of two and 2 m l < N, (5.17) can be simplified + m\ — c + (0) u(0) + ai^u(N) + T.l^{ Q2 r + 2iV-i }tl(l 1 2i£ 1 — a 2™ Series in expressions for c + (0), c(N — 1). and c(N) for filters from Figures 4.2(h) and 4.2(j) with m G [1. M — 1] can be, similar to (5.13) and (5.14), truncated according to the desired precision. For orders p greater than three, we implement B~ 1 {z) as a cascade of terms E{z) with different q's. Note that the output from block E(z) is always of the same type as the input to it. CHAPTER 6 APPLICATIONS In this chapter, we present a. comparison between multiscale spline derivative- based transform and traditional wavelet techniques on two sample applications. 6.1 Image Fusion Image fusion combines particular aspects of information from the same imag- ing modality or from distinct imaging modalities and can be used to improve the reliability of a particular computational vision task or to provide a human observer with a deeper insight about the nature of observed data. Whether it is combining different sensors or extending the dynamic range of a single sensor, the goal is to achieve more accurate inferences that can be achieved by a single sensor or a single sensor setting. The simplest method of fusing images is accomplished by computing their average. Such a technique does combine features from input images in the fused image, however, the contrast of the original features can be significantly reduced. Among more sophisticated methods, multiscale and multiresolution analyses have become particularly popular. Different pyramids [5, 41] and wavelet-based techniques [7, 20, 26. 32] have been applied to this problem. Fusion is performed in the transform space by computing local statistics across the decomposition scales. Typical size of neighborhood is between a single pixel and 5 x 5 area with some loss of contrast reported for the latter [5]. Some authors argue about of the particular criteria for fusion [5, 26], while others propose a set of different combination rules to perform image fusion in a variety of ways [7]. :»!) 60 Here, we are not trying to find the best criteria for fusion. We simply use a "maximum magnitude" rule (i.e., at each position and scale of the transform the corresponding transform coefficients are compared and those with greater magnitude included for reconstruction) and compare the results obtained by traditional wavelet techniques with the ones produced by the multiscale spline derivative- based trans- form. Figure 6. 1(a), (b) shows a pair of images with distinct (lower in Figure 6.1(a) and upper in Figure 6.1(b)) areas blurred. Figure 6.1(c) demonstrates the ideal result of image fusion obtained by manually cutting and pasting the upper part of Figure 6.1(a) and the lower part of Figure 6.1(b). Image fusion resulting from the use of discrete wavelet transforms is shown in Figure 6.1(d) for the orthogonal wavelet DAUB12 [7] and in Figure 6.1(e) for the linear phase biorthogonal wavelet Bior6.8. Figure 6.1(f) illustrates the result of image fusion with a multiscale spline derivative- based transform being employed. Note that the image fused using the multiscale spline derivative-based transform is virtually indistinguishable from the one obtained manually. Both orthogonal and biorthogonal wavelet transform exhibit artifacts. On the average, biorthogonal transform performs better than orthogonal due to the linear phase of wavelets and filters, but still causes serious artifacts. Basing the decision rule on an area rather than on a single pixel certainly reduces artifacts of traditional wavelet methods, however, it reduces sharpness as well. As noted by Burt and Kolczynski [5], "there is a tradeoff between sharpness and shift invariance." Here, we argue the point that the proper choice of transform can eliminate artifacts due to translation and rotation invariance altogether. Another possible source of artifacts is misregistration. Low-level fusion al- gorithms typically require input images that have already been properly aligned. Misregistration causes artifacts [7], and it can be recognized from our discussion 1,1 in Section 1.1 that they are more obvious for translation and rotation-noninvariant wavelet techniques than for the multiscale spline derivative-based transform. 6.2 Image Enhancement Similar to the previous section, we will point out problems associated with traditional wavelet analysis on an example. Below, we briefly present the problem of contrast enhancement in digital mammography. In mammography, early detection of breast cancer relies upon the ability to distinguish between malignant and benign mammographic features. Contrast en- hancement can make more obvious unseen or barely seen features of a mammogram without requiring additional radiation. The complexity of mammographic images and the subtlety of malignancies can present a challenge even to expert radiologists. In addition to dealing with barely visible mammographic features, human observers sometimes simply overlook abnormalities. Such mistakes affect the number of false- negative cases considerably. Computer processing of digital mammograms can assist radiologists to reach a correct diagnosis more consistently. A variety of computer based techniques have been reported in almost three decades of research. The advent of direct digital mammography devices has made digital image processing techniques more attractive for screening. Contrast enhancement in a wavelet transform framework can be achieved by applying a (typically nonlinear) function to wavelet coefficients and then reconstruct- ing an enhanced image with modified coefficients. We are not trying to find the best method for a specific type of images or noise here, rather we want to illustrate pos- sible sources of artifacts when orthogonal and biorthogonal wavelet transforms are used. Figure 6.2(a) shows an original mammographic image. Figure 6.2(b) depicts the result of enhancement using piecewise linear enhancement function [11, 23] to 62 (c) (d) Figure 6.1: (a) An image with lower part blurred, (b) An image with upper part blurred, (c) Fused image obtained by combining images from (a) and (b) man- ually. Results of fusion performed by (d) discrete orthogonal wavelet transform (DAUB12), (e) discrete biorthogonal wavelet transform (Bior6.8), and (f) multiscale spline derivative-based transform. 63 modify wavelet transform coefficients, and Figure 6.2(c) shows the result of enhance- ment using the multiscale spline derivative-based transform and the same enhance- ment function as Figure 6.2(b). In this figure, we. more than trying to point out better enhancement results in Figure 6.2(c). observe that Figure 6.2(b) shows fea- tures (artifacts) that are not present in the original image (they are particularly obvious in the area right from the mass in Figure 6.2(b)). Here, the goal of enhance- ment is making visible features that are present in the original image but obstructed by different structures within the same image. Generation of artificial features, while not being rare when nonlinearly processing coefficients of the orthogonal or biorthog- onal wavelet transform, is certainly not desirable for this particular application. For more details on different enhancement strategies in mammography, please refer to [11. 22, 23, 25]. Another problem related to enhancement and other kinds of processing of dig- ital mammograms (and not unique to mammography) is that the same mammogram could have been easily translated and/or rotated. From obvious reasons, it is not desirable that the result of processing can be significantly affected by positioning only — one more reason for a translation and rotation-invariant transform. Figure 6.3 shows outputs from steerable basis filters applied to the image from Figure 6.2(a) at scales {1.2.4}. By linearly combining these outputs analysis along arbitrary direction can be performed (cf. Chapter 4). This enables, as illustrated in Chapter 1. a rotation-invariant processing [21]. 64 Figure 6.2: (a) An original mammographic image containing a spicular mass. An enhanced image using (b) discrete biorthogonal wavelet transform (Bior6.8) and (c) multiscale spline derivative-based transform. (>.-> Figure 6.3: The magnitude of filter coefficients for G%'{u} T ,u<jy) at fc = {0, \t, \tc] (three levels of analysis shown). CHAPTER 7 CONCLUSION In this thesis we constructed a new transform that does not introduce artifacts due to translation and rotation invariance, which are inherent to traditional wavelet analysis. We extended the one-dimensional discrete dyadic wavelet transform to higher- order derivatives and even-order spline functions and developed an improved initial- ization procedure. Comparison to the originally employed initialization [28] showed significantly better performance of our procedure for finer scales of analysis. We developed several two-dimensional transforms. All of them were derived with the goal of eliminating artifacts due to lack of translation and rotation in- variance. We presented a two-dimensional discrete dyadic wavelet transform with a first-derivative wavelet as an extension of the one originally proposed in [28], a two- dimensional discrete dyadic wavelet transform with a second-derivative wavelet that can approximate Laplacian of a. Gaussian, and a steerable dyadic wavelet transform implemented in a near-perfect reconstruction filter bank which may be preferred when there is a need for orientation processing along equally spaced angles. We derived a multiscale spline derivative-based transform which uses x-y separable filters in a perfect reconstruction filter bank and enables fast translation and rotation-invariant directional analysis of images. We empirically showed the presence of artifacts in image fusion and enhance- ment applications when traditional wavelet methods were used. Here developed trans- forms did not exhibit similar artifacts, and, in case of fusion, yielded nice results without the use of postprocessing. 66 67 Future work will concentrate on a more thorough exploitation of the transform space for the existing applications, and potential new applications (e.g., denoising). In particular, we will try to use the causality property of a Gaussian kernel mentioned in Chapter I to trace features across scales and process them according to their multiscale behavior. REFERENCES 1] A. Aldroubi, M. Unser, and M. Eden, "Cardinal spline filters: Stability and convergence to the ideal sine interpolator," Signal Process., vol. 28, no. 2, pp. 127-138, 1992. 2] J. Babaud. A. P. Witkin, M. Baudin, and R. 0. Duda, "Uniqueness of the Gaussian kernel for scale-space filtering," IEEE Trans. Pattern Anal. Mach. Intel!., vol. 8. no. 1, pp. 26-33, 1986. 3] C. de Boor, A Practical Guide to Splines, Springer- Verlag, New York, NY, 1978. 4] P. J. Burt and E. H. Adelson, "The Laplacian pyramid as a compact image code," IEEE Trans. Commun., vol. 31, no. 4, pp. 532-540, 1983. 5] P. J. Burt and R. J. Kolczynski, "Enhanced image capture through fusion," in Proc. Int. Conf. Comput. Vision, Berlin, Germany, 1993, pp. 173-182. 6] J. Canny, "A computational approach to edge detection," IEEE Trans. Pattern Anal. Mach. IntelL, vol. 8. no. 6, pp. 679-698, 1986. 7] L. J. Chipman, T. M. Orr, and L. N. Graham, "Wavelets and image fusion," in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995, vol. 3, pp. 248-251. 8] C. K. Chui, An Introduction to Wavelets, Academic Press, San Diego, CA, 1992. 9] C. K. Chui, Ed., Wavelets: A Tutorial in Theory and Applications, Academic Press. San Diego, CA, 1992. [10] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. [\Y\ J. Fan and A. Laine, "Multiscale contrast enhancement and denoising in digital radiographs," in Wavelets in Medicine and Biology, A. Aldroubi and M. Unser, Eds.. CRC Press, Boca Raton, FL, 1996, pp. 163-189. [12] W. T. Freeman and E. H. Adelson, "The design and use of steerable filters," IEEE Trans. Pattern Anal. Mach. IntelL, vol. 13, no. 9, pp. 891-906, 1991. [13] H. Greenspan, S. Belongie, R. Goodman, P. Perona, S. Rakshit, and C. H. Anderson, "Overcomplete steerable pyramid filters and rotation invariance," in Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., Seattle, WA, 1994, pp. 222-228. [14] Y. Hel-Or and P. C. Teo, "Canonical decomposition of steerable functions," in Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., San Fran- cisco, CA, 1996, pp. 809-816. 69 [15] M. Holschneider, R. Kronland-Martinet, J. Morlet, and Ph. Tchamitchian, "A real-time algorithm for signal analysis with the help of the wavelet transform," in Wavelets, Time-Frequency Methods and Phase Space, J. M. Combes, A. Gross- mann, and Ph. Tchamitchian, Eds., Springer- Verlag, Berlin, Germany, 1989, pp. 286-297. [16] A. J. Jerri, "The Shannon sampling theorem — its various extensions and appli- cations: A tutorial review," Proc. IEEE, vol. 65, no. 11, pp. 1565-1596, 1977. [17] A. Karasaridis and E. Simoncelli, "A filter design technique for steerable pyramid image transforms," in Proc. IEEE Int. Conf Acoust. Speech Signal Process., Atlanta, GA, 1996, vol. 4, pp. 2389-2392. [18] I. Koren and A. Laine, "A discrete dyadic wavelet transform for multidimensional feature analysis," in Time-Frequency and Wavelet Transforms in Biomedical Engineering, M. Akay, Ed.. IEEE Press, New York, NY, 1997. [19] I. Koren, A. F. Laine, J. Fan, and F. J. Taylor, "Edge detection in echocardio- graphic image sequences by 3-d multiscale analysis," in Proc. IEEE Int. Conf. Mage Process., Austin, TX, 1994, vol. 1, pp. 288-292. [20] I. Koren. A. Laine. and F. Taylor, "Image fusion using steerable dyadic wavelet transform," in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995, vol. 3, pp. 232-235. 21}) I. Koren, A. Laine. F. Taylor, and M. Lewis, "Interactive wavelet processing and techniques applied to digital mammography," in Proc. IEEE Int. Conf. Acoust. Speech Signal Process., Atlanta, GA, 1996, vol, 3, pp. 1415-1418. ] A. Laine, J. Fan, and S. Schuler, "A framework for contrast enhancement by dyadic wavelet analysis." in Digital Mammography, A. G. Gale, S. M. Astley, D. R. Dance, and A. Y. Cairns, Eds., Elsevier, Amsterdam, The Netherlands, 1994, pp. 91-100. A. Laine, J. Fan, and W. Yang, "Wavelets for contrast enhancement of digital mammography," IEEE Eng. Med. Biol. Mag., vol. 14, no. 5, pp. 536-550, 1995. [24] A. Laine, I. Koren, W. Yang, and F. Taylor, "A steerable dyadic wavelet transform and interval wavelets for enhancement of digital mammography," in Wavelet Applications II. Proc. SPIE, H. H. Szu, Ed., Orlando, FL, 1995, vol. 2/91, pp. 736-749. [25] A. F. Laine, S. Schuler, J. Fan, and W. Huda, "Mammographic feature en- hancement bv multiscale analysis," IEEE Trans. Med. Imaging, vol. 13. no. 4, pp. 725-740, 1994. [26] H. Li, B. S. Manjunath, and S. K. Mitra, "Multi-sensor image fusion using the wavelet transform," in Proc. IEEE Int. Conf. Image Process., Austin, TX, 1994, vol. 1, pp. 51-55. [27] S. Mallat and W. L. Hwang, "Singularity detection and processing with wavelets," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 617-643, 1992. [28] S. Mallat and S. Zhong, "Characterization of signals from multiscale edges," IEEE Trans. Pattern Anal. Mach. Inteli, vol. 14, no. 7, pp. 710-732, 1992. 71) [29] D. Marr and E. Hildreth, "Theory of edge detection," Proc. R. Soc. London Ser. B, vol. 207, no. 1167. pp. 187-217, 1980. [30] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, Prentice- Hall, Englewood Cliffs, NJ, 1989. [31] P. Perona. "Deformable kernels for early vision," IEEE Trans. Pattern Anal. Mach. Intel!., vol. 17, no. 5, pp. 488-499, 1995. [32] T. Ranchin, L. Wald, and M. Mangolini, "Efficient data fusion using wavelet transform: the case of SPOT satellite images," in Matin nmtical Imaging: Wavelet Applications in Signal and Image Processing. Proc. SPIE. A. F. Laine, Ed., San Diego, CA, 1993, vol. 2034, pp. 171-178. [33] 0. Rioul and P. Duhamel, "Fast algorithms for discrete and continuous wavelet transforms," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 569-586, 1992. [34] I. J. Schoenberg, "Contributions to the problem of approximation of equidistant data by analytic functions," Q. Appl. Math., vol. 4. no. 1. pp. 45-99, 112-141, 1946. ' [35] I. J. Schoenberg. "Cardinal interpolation and spline functions." ./. Approx. Theory, vol. 2, no. 2, pp. 167-206, 1969. [36] I. J. Schoenberg, "Notes on spline functions III: On the convergence of the interpolating cardinal splines as their degree tends to infinity," Isr. J. Math.. vol. 16. no. 1, pp. 87-93, 1973. [37] C. E. Shannon, "Communication in the presence of noise," Proc. IRE. vol. 37, no. 1, pp. 10-21, 1949. [38] E. P. Simoncelli and W. T. Freeman, "The steerable pyramid: A flexible archi- tecture for multi-scale derivative computation," in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995, vol. 3, pp. 444-447. [39] E. P. Simoncelli, W. T. Freeman, E. H. Adelson, and D. J. Heeger, "Shiftable multiscale transforms," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 587-607, 1992. [40] M. J. T. Smith and T. P. Barnwell, III, "Exact reconstruction techniques for tree-structured subband coders," IEEE Trans. Acoust. Speech Signal Process., vol. 34, no. 3, pp. 434-441, 1986. [41] A. Toet, "Image fusion by a ratio of low-pass pyramid." Pattern Recognit. Lett., vol. 9. no. 4, pp. 245-253, 1989. [42] M. Unser and A. Aldroubi, "A general sampling theory for nonideal acquisition devices," IEEE Trans. Signal Process., vol. 42, no. 11, pp. 2915-2925, 1994. [43] M. Unser, A. Aldroubi, and M. Eden, "Fast B-spline transforms for continuous image representation and interpolation," IEEE Trans. Pattern Anal. Mach. Intel!., vol. 13. no. 3, pp. 277-285, 1991. [44] M. Unser, A. Aldroubi. and M. Eden, "On the asymptotic convergence of B- spline wavelets to Gabor functions," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 864-872, 1992. 71 [45] M. Unser, A. Aldroubi, and M. Eden, "Polynomial spline signal approximations: Filter design and asymptotic equivalence with Shannon's sampling theorem," IEEE Trans. Inf. Theory, vol. 38, no. 1, pp. 95-103, 1992. [46] M. Unser, A. Aldroubi, and S. J. Schiff, "Fast implementation of the continuous wavelet transform with integer scales," IEEE Trans. Signal Process., vol. 42, no. 12, pp. 3519-3523, 1994. [47] M. Vrhel, C. Lee, and M. Unser, "Fast computation of the continuous wavelet transform through oblique projections," in Proc. IEEE Int. Conf. Acoust. Speech Signal Process., Atlanta, GA, 1996, vol. 3, pp. 1459-1462. [48] A. Witkin. "Scale space filtering," in Proc. Int. Joint Conf. Artif. IntelL, Karlsruhe, Germany, 1983, pp. 1019-1022. BIOGRAPHICAL SKETCH Iztok Koren earned the B.S. and M.S. degrees in electrical engineering from the University of Ljubljana, Slovenia, in 1987 and 1991, respectively. He will receive the Ph.D. degree in electrical and computer engineering from the University of Florida, Gainesville, in 1996. His research interests include image processing, digital signal processing, and wavelets. He is a member of Eta Kappa Nu, Tau Beta Pi, and the IEEE. 7-2 I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. h/4/ /f4f/fr^ ,Fred J/Taylor/, Chairman Profess'or of/Electrical and Computer Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. 6^/iU. yj , -<^£Ui> Andrew F. Laine , Cochairman Associate Professor of Computer and Information Science and Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and-^-fylly adequate, in scope and quality, as a dissertation for the degree of Doctor / Jose C/ olr of Ele^rical and Computer ngmeering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of DocJ*^ of Philosophy. 'John M. M. Anderson Assistant Professor of Electrical and Computer Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doyetor of Philosophy. ■pai 6/lA ?mit Sigmon Associate Professor of Mathematics This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1996 V^ Winfred M. Phillips Dean, College of Engineering Karen A. Holbrook Dean, Graduate School L LD 1780 19S£ .Km UNIVERSITY OF FLORIDA 3 1262 08555 0696